//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

// Generate bfield 3D density map (voxel data) in situs format

#include "ephi.hpp"

#define RADIUS 0.15
#define WR 0.035
#define SPACING 0.01

#define GRIDN 250
#define CONTACT ((float)9e99)

static prec_t map_range;
static prec_t voxel_size;

struct task : Task
{
  task(Statics& statics, size_t z, float *wp) : statics(statics), z(z), wp(wp), maxf(0) { }
  Statics& statics;
  size_t z;
  float *wp;
  float maxf;

  void execute ()
  {
    float *o = wp;
    for ( size_t yi = 0; yi < GRIDN; yi++ )
      {
        for ( size_t xi = 0; xi < GRIDN; xi++ )
          {
            vect3d pos(-map_range + xi * voxel_size, -map_range + yi * voxel_size, -map_range + z * voxel_size);
            float mag;
            if ( statics.inContact (pos) )
              {
                mag = CONTACT;
              }
            else
              {
                vect3d bf, ef;
                statics.getFields (pos, bf, ef);
                mag = prec2double(bf.length());
                if ( maxf < mag ) maxf = mag;
              }
            *o++ = mag;
          }
      }
  }
};

int main(int argc, char *argv[])
{
  Statics statics;
#if 1
  make_polywell_cube (statics, RADIUS, WR, SPACING, 1, 1e-7);
#else
  CoilData cd;
  cd.load ("out.cgen");
  make_polywell_cube (statics, RADIUS, WR, SPACING, 1, 1e-7, CoilDataFactory<StaticCIL>(cd));
#endif

  prec_t rr = RADIUS + (WR * 2 + SPACING) / PREC_SQRT2;
  map_range = rr + 2 * WR;
  voxel_size = 2 * map_range / GRIDN;

  FILE *fp = fopen ("mapgen.sit", "w");
  // voxel size (float), origin xyz (float), grid size xyz (int)
  fprintf (fp, "%e %e %e %e %d %d %d\r\n", prec2double(voxel_size), 0.0, 0.0, 0.0, GRIDN, GRIDN, GRIDN);

  size_t sc = GRIDN * GRIDN * GRIDN;
  float *data = new float [sc];
  TaskManager tmgr;
  for ( size_t zi = 0; zi < GRIDN; zi++ )
      tmgr.addTask (new task (statics, zi, data + zi * GRIDN * GRIDN));
  float maxf = 0;
  task *t;
  while ( (t = (task *)tmgr.run()) )
    {
      if ( maxf < t->maxf )
          maxf = t->maxf;
      delete t;
    }
  //exit(0);

  bool had_nl = true;
  for ( size_t i = 0; i < sc; i++ )
    {
      float v = data[i];
      if ( v == CONTACT )
          v = maxf;

      if ( had_nl )
        {
          fprintf (fp, "%.3e", v);
          had_nl = false;
        }
      else
        {
          fprintf (fp, " %.3e", v);
        }

      if ( i && !(i % 128) )
        {
          fprintf (fp, "\r\n");
          had_nl = true;
        }
    }
  fclose (fp);

  delete [] data;

  return 0;
}

